Emergent self-organized complex network topology out of stability constraints 
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Although most networks in nature exhibit complex topologies the origins of such complexity 
remain unclear. We propose a general evolutionary mechanism based on global stability. This 
mechanism is incorporated into a model of a growing network of interacting agents in which each 
new agent's membership in the network is determined by the agent's effect on the network's global 
stability. It is shown that out of this stability constraint, complex topological properties emerge in 
a self organized manner, offering an explanation for their observed ubiquity in biological networks. 
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Complex networks of interacting agents are ubiquitous, 
in a wide range of scales, from the microscopic level of 
genetic, metabolic and proteins networks to the macro- 
scopic human level of the Internet [l|, 0|. All of them 
exhibit high clustering and relatively short path length 
compared with random networks. In addition, they fre- 
quently show a nonhomogeneous structure, characterized 
by a degree distribution (the probability of a node to be 
connected to k other ones) with a broad tail Pjk)^ k"'^ 
for large values of k, with exponents 7 < 3 [l|, d, S| • Sev- 
eral mechanisms have been proposed to give rise to this 
kind of topologies [1, 4]. These mechanisms have success- 
fully explained the origin of complexity in some networks, 
but it is recognized that another, equally large, number 
of cases can not be accounted for either class of mod- 
els. In particular, growing biological networks involve 
the coupHng of at least two dynamical processes. The 
first one concerns the addition of new nodes, attached 
either during a slow evolutionary (i.e., species lifetime) 
or a relatively faster developmental (i.e., organism life 
time) process. A second one is the node dynamics which 
affects and in turn is affected by the growing processes. 
It is reasonable to expect that the network topologies we 
finally witness could have emerged out of these coupled 
processes. This Letter is dedicated to discuss a simple 
model of this problem, showing that complex networks 
do emerge under general realistic constraints. It needs 
to be noted from the outset that the aim of this letter 
is not to describe an arbitrary algorithm, but to identify 
a dynamical process able to be implemented by natural 
systems. 

Before introducing the model, and to fix ideas, let us 
dwell on some concrete general examples. First con- 
sider a food web, which is constructed through commu- 
nity assembly rules, strongly influenced by the underly- 
ing dyriamics of species and speciflc interactions among 
themes, Another example could be neuronal networks, 
where the addition of hundreds of thousands of new neu- 
rons is followed by a dynamical process in which neuronal 



dynamics and connectivity are interrelated in a way not 
fully understood. Yet a third example at another scale, 
could be imagined in the context of social networks, in 
which novice members can be accepted or rejected based 
on their individual contribution to a global interest, flt- 
ness, performance or profit. In the three examples it is 
relatively easy to visualize the two processes mentioned 
above. The consequence of adding a new member with a 
given connectivity affecting a global in / stability, is repre- 
sented in these examples by the aboundance/lack of food, 
the neuronal welfare/death or the profits' up/down, re- 
spectively. Notice that each new member may not only 
result in its own addition/rejection to the system, but it 
can also promote avalanches of extinctions amongst exist- 
ing members, an effect we found that strongly infiuences 
the network's topology. 

Let us consider a system of n interactive agents, 
whose dynamics is given by a set of differential equa- 
tions dx/dt — F{x), where x is an n-component vector 
describing the relevant state variables of each agent and 
F is an arbitrary non-linear function. One could imagine 
that X in different systems may represent concentrations 
of some hormones, or the average density populations 
in a food web, or the concentration of a chemicals in a 
biochemical network, or the activity of genes in a gene 
regulation net, etc. We assume that a given agent i in- 
teracts only with a limited set of fc^ < n other agents; 
thus Fi depends only on the variables belonging to that 
set. This defines the interaction network, as was done 
previously 0. 

We will assume that there are two time scales in the dy- 
namics. On the long time scale (much larger than the ob- 
servation time) the system is subjected to an external flux 
(migration, mutation, etc.) of new agents that interact 
with some of the previous ones and can be incorporated 
into the system or not, so n (and the whole set of differ- 
ential equations) can change. On short time scales we as- 
sume that n is constant and the dynamics already led the 
system to a particular stable stationary state x* deflned 
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by F{x*) = Q. The stability of that solution is deter- 
mined by the eigenvalue with maximum real part of the 
Jacobian matrix aij = ( ) . Therefore a new agent 

will be incorporated to the network if its inclusion result 
in a new stable fixed point, that is, if the values of the 
interaction matrix j are such that the eigenvalue with 
maximum real part A of the enlarged Jacobian matrix 
is negative (A < 0). Assuming that isolated agents will 
reach stable states by themselves after certain character- 
istic relaxation time, the diagonal elements of the matrix 
Gi^i are negative and given unity value to further simplify 
the treatmentj3|. The interaction values, (i.e., the non- 
diagonal matrix elements aij) will take random values 
(both positive and negative) taken from some statisti- 
cal distribution. In this way we have an unbounded en- 
semble of systems [3| characterized by a "growing through 
stability" history. Randomness would be self-generated 
through the addition of new agents processes. Each spe- 
cific set of matrix elements after addition defines a par- 
ticular dynamical system and the subsequent analysis for 
time scales between successive migrations is purely deter- 
ministic. 

These ideas are implemented in a numerical model as 
follows: At every step the network can either grow or 
shrink. In each step an attempt is made to add a new 
node to the existing network, starting from a single agent 
{n = 1). Based on the stability criteria discussed, the at- 
tempt can be successful or not. If successful, the agent 
is accepted, so the existing n x n matrix grows its size 
by one column and one row. Otherwise the novate agent 
will have a probability to be deleted together with some 
other nodes as further explained below. More specifically, 
suppose that we have an already created network with n 
nodes, such that the n x n associated interaction matrix 
Gij is stable. Then, for the attachment of the n + 1th 
node we first choose its degree kn+i randomly between 
1 and n with equal probability. Then the new agent in- 
teraction with the existing network member i is chosen 
such that non-diagonal matrix elements {ui^n+ii <in+i,i) 
(i — 1, . . . , n) are zero with probability 1 — fc„+i/n and 
different from zero with probability fc„+i/n; to each non- 
zero matrix element we assign a different real random 
value uniformly distributed in [—5,6]. b determines the 
interaction range variability and it is one of the two pa- 
rameters of the model. Then, we calculate numerically A 
for the resulting (ri,-|-l)x(n-|-l) matrix. K A < the new 
node is accepted. If A > it means that the introduction 
of the new node destabiHzed the entire system and we 
will impose that, either the new agent is eliminated or it 
remains but produces the extinction of a certain number 
of previous existing agents. In order to further simplify 
the numerical treatment, we will allow up to <? < kn+i ex- 
tinctions, taken from the set of kn+i nodes connected to 
the new onef^; q is the other parameter of the model. To 
choose which nodes are to be eliminated, we first select 
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FIG. 1: (Color on-line) Degree distribution P(k) for b = 2, 
q = 3 and different values of Umax; the dashed lines corre- 
spond to a power law P{k) ~ k ' with 7 « 2.4. The inset 
shows 7 as a function of b for different values of q. 



one with equal probability in the set of fc„+i and remove 
it. If the resulting n x n matrix is stable, we start a 
new trial; otherwise, another node among the remaining 
kn+i — 1 is chosen and removed, repeating the previous 
procedure. If after q removals the matrix remains unsta- 
ble, the new node is removed, we return to the original 
n X n matrix and start a new trial[lo|. 

First we calculated the average connectivity C(n), de- 
fined as the fraction of non-diagonal matrix elements dif- 
ferent from zero, averaged over different runs. We found 
that C{n) ~ (see Supplementary information) for 

large values of n, where the exponent e depends on b and 
q, taking values < e < 1. Such behavior is characteristic 
of food webs[_ll|| and it has been interpreted in terms of 
self-organized criticality concepts[l3|; the present results 
suggest that this is a general behavior in stability-driven 
self organized systems. 

Next we calculated the degree distribution P{k) of the 
network with n = Umax for different values of b and q. 
The typical behavior of P{k) is illustrated in Fig[T] for 
6 = 2, g = 3 and different values of rimax- We see the 
emergence of a fat tail P{k) ^ k^"' for large values of n, 
with an exponent 7, independently of the network size 
(this figure also shows that the drop in the tail of the 
distribution is a finite size effect). Notice that this rela- 
tively small range of the broad tail is what more often is 
seen in real networks. The qualitative behavior of P{k) 
for other values of b and q is the same (see Supplemen- 
tary information) . The inset of Fig[T] shows the value of 
the exponent 7 as a function of b for different values of 
q. We see that 7 presents a minimum around 6 = 2 for 
all values of q; as q increases the exponent decreases and 
for large enough values of q we obtain a non-trivial value 
of 7 < 3 for a broad range of values of 6. 

To exclude the possibility that the observed network 
topology is trivially associated with a hidden preferen- 
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FIG. 2: (Color on-line) Relative attachment probability 
Il{k) / P{k) for different values of Umax, compared with the 
corresponding results for a BA network of the same size. 
Dashed line corresponds to a linear behavior Il{k)/ P(k) ~ k. 



tial attachment (PA) process, we computed the attach- 
ment probability Tl(k), defined as the probabiHty that a 
new node connects with an aheady existing node with 
degree k. Assuming that the average degree (ki) ^ n, 
the attachment probabiUty can be expressed as n(fc) = 
Hi, where 11^ is the probability that the new node 
connects to the already existing node i, Uk « nP{k) is 
the number of nodes with degree k and the sum runs over 
all sites i with degree ki — k. If stability selection would 
favor some kind of PA mechanism, (i.e., if new nodes 
are attached with larger probability to nodes highly con- 
nected) we should expect 11^ 
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FIG. 3: (Color on-line) (a) Networks' average clustering co- 
efficient Cc(n) and L{n) fov b — 2 and different values of q 
as functions of network size. The dashed line is a guide to 
the eye corresponding to Cc{n) ^ n""'^^. Inset shows the 
same data plotted against each other, (b) Cc(n) and L{n) as 
functions of network size computed from a single network re- 
alization. Data are samples taken every fifty trials, regardless 
of the resulting stability. Notice how fluctuations increase as 
the network grows. Inset shows the same data plotted against 
each other (full circles), in addition to the data computed from 
a random network with equal size and density of connections 
(open circles). 



In Fig. [2] the relative attachment probability 
Ii{k) / P{k) in the present model for a fixed network size 
n and diflerent values of b is compared with the corre- 
sponding results for a network of Jhe same size obtained 
with the Barabasi-Albert (BA) [li] algorithm with con- 
nectivity C{n). This quantity shows the expected be- 
havior n(A:)/P(fc) ^ k for large values of fc, consistently 
with Eq.|(T|). In the present model Ii{k)/P{k) remains 
almost constant for a wide range of values of k (includ- 
ing a range of values for which the power law behavior 
of P(fc) has already established), but displays an increas- 
ing trend consistent with Eq.(IT]) for large values of k. In 
other words, in the present model at variance with the 
BA model, as the network grows, the assembly mecha- 
nisms selected by stability shows a crossover between two 
regimes: one dominated by PA and the other not. 

Considering that biological systems are probably never 
in a completely stable situation, we relaxed the condition 
of stability A < and look at networks growing by allow- 



ing A to take small positive values so that the characteris- 
tic time to leave an unstable fixed point t — A^^ » 1. By 
accepting nodes as long as A < A the calculation of P{k) 
for different values of A (positive and negative) showed 
similar qualitative behavior, with small variations of the 
7 exponent (See Supplementary Information). 

Next we calculated the average path length L between 
two nodes and the average cluster coefficient Cc for the 
networks obtained by the present algorithm as a function 
of the network size n. L is defined as the minimum num- 
ber of links needed to connect any pair of nodes in the 
network and Cc is defined as the fraction of connections 
between topological neighbors of any site[ll|. In Fig. [3] 
we show the typical behavior of L{n) and Cc{n). We see 
that Cc{n) ^ n^^-"^^ and L{n) ~ A Inn+C. Such scaling 
behavior is the same one observed in the BA modelQ. 

As shown in Fig. [3] larger networks becomes less clus- 
tered and have longer minimum path on the average. L 
and Cc are inversely related as it can be seen in the single 
run plotted in panel b of Fig. [3l The data correspond to 
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FIG. 4: (Color on-line) Degree distribution P{k) for fe = 2, 
q = 3, Umax = 100 and different values of R; the straight lines 
correspond to a power law fits. 



values computed every fifty trials, whether or not the at- 
tempt to add a node was successful or not at that trial. In 
a sense, this is how a natural network would look Hke to 
an observer if one could take snapshots in time. Clearly 
both quantities fluctuate in opposite directions, as fur- 
ther shown in the inset where the data corresponding to 
a randomly shuffled network is also plotted for compari- 
son. The behavior of Cc and L is Hnked with the selection 
dynamics ruling which node is accepted or rejected. The 
stability constraint favors the nodes with few links, since 
they modify the matrix aij stability much less than new 
nodes with many links (of course this is reflected in the 
P{k) density). Thus, most frequently the network grows 
at the expense of adding nodes with one or few links, 
producing an increase of L and a decreases of Cc. Most 
of the times, nodes with many links destabilize the net- 
work and are rejected, but when one is flnally accepted, 
a large decrease in L together with an increase in Cc is 
observed. This sudden change is the signature of a new 
network hub, as seen in the example denoted with an ar- 
row in Fig. 4b. We also verifled that those fluctuations 
lead to a slow diffusive-like growth of the network size 
n{t) ~ t^/^ (not shown), where the time is measured in 
number of trials. 

Finally, to consider the effect of local stability selection 
pressure, we modifled the algorithm as follows. Once the 
new candidate node and its 1st nearest neighbors (nn) 
are chosen, we analyze the stability of the subnetwork 
composed by neighbors up to a range R {R = 1: 1st nn, 
R = 2: 2nd nn, etc.). In Fig [4] we show P{k) at n — 
n-max for different values of R. We see that the fat tail 
P(fc) ~ k~'^ appears as long as i? > 3, which coincides 
with the value of L for the corresponding net size (see Fig. 
3) , suggesting a correlation between stability and the self 
organized emergence of small world topology. Notice also 
that considering local stability allows a larger variability 



in the value of 7 (7 « 0.9 for i? = 3), although 7 quickly 
converges to the global stability result (for i? > 4 both 
results are almost indistinguishable). 

The robustness of complex networks against error and 
attack has already been investigated[l3| considering the 
effects of nodes or links' deletion. The present results 
shows that the consequences of perturbing a single node 
may depend also on stability, a topic that deserves further 
clariflcation. Indeed, closely related results on Boolean 
networks dynamics supports already the generality of this 
approach [U]. 

Summarizing, the analysis of a simple model shows 
that complex topology can appear in networks as an 
emergent property driven by a stability selection pres- 
sure during the growing process. This suggests yet an- 
other explanation for the ubiquity of complex topology 
observed in different networks in nature. 
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